Introduction

Forests are a defining feature of the US North. This region is the most heavily forested region in the US, and these forests are an essential provider of employment and recreation opportunities and other ecosystem services for the region’s people, and of habitat for the region’s wildlife (Shividenko et al. 2005, Shifley et al. 2012). In addition, the forests of the Upper Midwest in particular have been a net carbon sink for most of the past century, as fast-growing early successional species established and thrived following near-complete deforestation around the end of the 19th Century (Birdsey et al. 2006, Williams et al. 2012). However, the fate of these forests is highly uncertain. For one, many Upper Midwest forests are undergoing a shift in community composition from early-successional to mid-successional species, and the ecological and biogeochemical consequences of this transition are not well-established (Gough et al. 2016). In addition, these forests face a multitude of direct and indirect pressures from both biotic and abiotic sources (Shifley and Moser 2016), including longer and more severe droughts (Gustafson et al. 2016, Liénard et al. 2016, Swanston et al. 2017), non-native insects and pathogens (Lovett et al. 2016), and more. Collectively, these processes have culminated in forests that are different from their pre-settlement climax counterparts (Wolter and White 2002, Stearns and Likens 2002, Thompson et al. 2013), and our ability to make predictions about such non-analog environments based on the past is limited.

More reliable predictions can in theory be obtained by using dynamic vegetation models that explicitly represent processes involved in forest growth and mortality. Vegetation models fall broadly along the following spectrum of complexity, which is a direct consequence of trade-offs between realism and process fidelity on one hand and tractability, computational demand, and data requirements on the other (Hawkes 2000). On one side are relatively simple “big leaf” models—such as PNET, SiBCASA, and Biome-BGC—in which the vegetation at a particular location consists essentially of a single large “plant” whose characteristics are the (weighted) average of all the vegetation at that site. A common application of these models was simulating the land surface boundary condition for atmospheric general circulation models, though similar models have been (and continue to be) applied successfully to simulate biogeochemical processes in many ecosystem science contexts. On the other end of the spectrum are “individual-based” models (a.k.a. “gap models”), which explicitly simulate multiple indivduals competing for resources at a single site (Shugart et al. 2015). Examples of such models are LANDIS (Mladenoff 2004, Scheller et al. 2007) and UVAFME. Because these models can explicitly represent inter-specific differences in plant productivity, resource allocation, and stress tolerance (among others), as well as the competition that emerges out of these differences, they may be better able to represent changes in ecosystem-scale processes, especially in no-analog conditions (Purves and Pacala 2008). Between these two extremes are models that use approximations attempt to capture the emergent biogeochemical and ecological outcomes of interactions between individual plants without the need for explicitly simulating each individual (Purves and Pacala 2008). One specific example of this approach is the Ecosystem Demography model (ED2; Moorcroft et al. 2001, Medvigy and Moorcroft 2011, Longo et al. 2019b), which is the focus of this study.

Regardless of vegetation model complexity, their projections are inherently uncertain. This uncertainty comes from many different sources, which can broadly be classified into the following categories: Driver uncertainty refers to uncertainty in data about processes not represented by the model (e.g. weather and climate for vegetation models). Initial condition uncertainty arises from the fact that models have to start somewhere, and the exact conditions at the place and time simulations begin are frequently unknown. Process or structural uncertainty arises because vegetation models necessarily only represent a subset of all processes involved in plant biology, and that the process that are included can usually be represented in multiple different ways (i.e. which processes are included, and how they are represented). Finally, parameter uncertainty arises because of natural variability and imperfect calibration of the above process representations. The relative importance of these sources of uncertainty in ecological forecasts is a key question because it is directly related to future research priorities. For example, in the atmospheric science community, the insight that meteorological forecasts were most sensitive to initial condition uncertainty (Lorenz REF) led to a multi-decadal research agenda aimed at constraining the initial state with improved observations, which has directly contributed to a steady and persistent improvement in meteorological forecasts over the last several decades.

Many ecosystem modeling studies have looked at these categories of uncertainty independently. Perhaps the most work has been done on parameter uncertainty in ecosystem models. Dietze et al. (2014) investigated which parameters were most important to ED2 predictions of C sequestration across a diverse range of sites in North America. The found that the parameters contributing the most to overall predictive uncertainty were those related to growth respiration, mortality, stomatal conductance, and water uptake, and that parameter uncertainty varied with biome. More recently, Raczka et al. (2018) evaluated parameter uncertainty over 100 year time scale at the Willow Creek Ameriflux site, and found that the most important parameters to model uncertainty were quantum efficiency of photosynthesis, leaf respiration, and soil-plant water transfer. Another recent study assessed parameter uncertainty in ED2 as part of its overall goal of parameterizing the model for sagebrush (Pandit et al. 2018). They found that the most important parameters were specific leaf area, the maximum rate of carbon fixation (Vcmax), slope of stomatal conductance, the fine-root to leaf carbon ratio, and the fine root turnover rate. Another study looked specifically at contributions to uncertainty from parameters related to canopy radiative transfer, and found that parameters related to both canopy structure and leaf optical properties had a large impact on predictions of net primary productivity (Viskari et al. 2019). {Driver uncertainty – a few examples?} Other studies have investigated structural uncertainty.. Numerous vegetation model intercomparison studies have demonstrated that different models produce significantly different projections of overall land carbon sequestration (Friedlingstein et al. 2006, 2014), response to CO2 fertilization (Zaehle et al. 2014, Walker et al. 2015, Medlyn et al. 2015), and soil C sequestration (Sulman et al. 2018), among others. In particular, Lovenduski and Bonan (2017) found that model structural uncertainty was the primary source of uncertainty in Earth system model projections of the land C sink. This structural uncertainty has been attributed to differences in model representations of key processes, including canopy radiative transfer (Fisher et al. 2017), soil biogeochemistry (Wieder et al. 2017, Sulman et al. 2018), stomatal conductance (Franks et al. 2018), and photosynthesis (Rogers et al. 2016). However, the extent to which specific processes contribute to model uncertainty is difficult to evaluate from model intercomparisons because different models are different from each other in too many different ways, meaning that structural effects are confounded with other uncertainty sources. Moreover, comparative analysis of contributions of different types of uncertainty to model projections are rare in the ecosystem modeling literature.

Previous work found that an earlier version of ED had difficulty sustaining multiple plant functional types in long simulations (Fisher et al. 2010). One of the primary drivers of successional change in temperate forests is competition for light , and differences in plants’ strategies for light acquisition and utilization are a key driver of niche differentiation in these ecosystems . To this end, we focused our structural uncertainty analysis on three processes related to competition for light: Representation of horizontal competition; representation of canopy radiative transfer (especially of diffuse radiation); and trait plasticity with light level.

This study focuses on the interaction between parametric and structural uncertainty in the ED2 model. Our study is organized around the following guiding questions: (1) Which processes related to light utilization are most important to accurately modeling community succession and C cycle dynamics in temperate forests of the Upper Midwest? (2) What is the cost of considering these processes, in terms of additional parametric uncertainty? (3) What are the relative contributions of parametric uncertainty (data limitation) vs. structural uncertainty (theoretical limitation?)? To answer these questions, we ran ED2 ensemble simulations with a factorial combination of submodels related to radiative transfer formulation (two-stream vs. multiple scatter), horizontal competition (finite canopy radius vs. complete shading), and trait plasticity (whether or not SLA and Vcmax vary with light level) We then used sensitivity and variance decomposition analyses to evaluate the contribution of parameter uncertainty for each model configuration, exploring the ecological implications of the results.

Methods

Site description

We performed this study at the University of Michigan Biological Station (UMBS; Ameriflux site US-UMd, 45.5625\(^{\circ}\), -84.6975\(^{\circ}\)), located in Northern Lower Michigan, USA. The area surrounding the research station is 87% well-drained upland forest and 13% wetland (Bergen and Dronova 2007); the focus of this study is on the former. The landscape geography of the UMBS upland forest is 20.4% moraine, 37.8% high outwash plain, 31.3% low outwash plain, 5.7% lake-margin terrace, 3.6% ice-contact, and 1.2% lowland glacial lake (Bergen and Dronova 2007). Most of the UMBS upland forest canopy is dominated by temperate deciduous early-successional species, most importantly Populus grandidentata (bigtooth aspen) and, to a lesser extent, Betula papyrifera (paper birch), with Acer saccharum (sugar maple), Acer rubrum (red maple), Fagus grandifolia (American beech), Tilia americana (basswood), Betula alleghaniensis (yellow birch), Fraxinus americana (white ash), Tsuga canadensis (eastern hemlock), Quercus rubra (northern red oak), Pinus strobus (white pine), and Pinus resinosa (red pine) existing in various fractions in the understory (and, in patches, in the canopy). This composition is a legacy of the site’s disturbance history: The site was intensively logged in the late 1800s and early 1900s, and experienced regularly recurring fires until the mid-1920s, at which point a regime of active fire suppression started that has persisted to the present day. As a result, the average stand age in 2013 was 95 years. The majority of the forest that is aspen-dominated is undergoing succession to “northern hardwood” (maple, beech, basswood, birch, ash, hemlock), “upland conifer” (red and white pine), or “northern red-oak” ecotypes (Bergen and Dronova 2007).

Ecosystem Demography Model (version 2.2)

For this study, we used the Ecosystem Demography Model, version 2.2 (ED-2.2). A full description of the default configuration of the model is provided by Longo et al. (2019b). Briefly, ED-2.2 solves the energy, water, and carbon cycles separately for each of multiple “cohorts” of trees of similar composition, size, and age sharing a micro-environment and disturbance history (Moorcroft et al. 2001, Medvigy et al. 2009, Longo et al. 2019b). Besides the cohort-based approximation of gap dynamics, a distinctive feature of ED-2.2 is that the equations for energy, water, and carbon fluxes are defined in terms of total energy, water, and carbon, which ensures excellent conservation of mass in long-running (multidecadal) simulations (Longo et al. 2019b). Various versions of the ED model have been validated in boreal (Medvigy and Moorcroft 2011), temperate (Medvigy et al. 2009, Dietze et al. 2014, Raczka et al. 2018), and tropical biomes (Moorcroft et al. 2001, Longo et al. 2019a), and have been applied in, among others, paleoecological studies (Rollinson et al. 2017), free-air CO2 enrichment studies (De Kauwe et al. 2013, Miller et al. 2015), estimation of potential carbon stocks (Hurtt et al. 2004), and analyses of regional vegetation-climate feedbacks (Swann et al. 2015).

The official ED2 source code is available at https://github.com/EDmodel/ED2, and the exact version used in this study (which includes minor revisions to accommodate the requirements of this study) can be obtained at https://github.com/ashiklom/ED2/tree/b048950971e91699b78fc566df133caf977fda32.

Submodels

In this study, we ran a factorial combination of the following ED-2.2 configurations: (1) two-stream vs. multiple-scatter canopy radiative transfer models; (2) infinite vs. finite crown area (a.k.a. complete vs. partial shading); and (3) static vs. light-plastic traits.

Both canopy radiative transfer models in ED-2.2 resolve the full vertical radiation profile within a patch as a function of canopy structure (leaf and wood area indices, crown area, leaf angle distribution) and incident solar radiation, following the definitions in CLM 4.5 (Oleson et al. 2013). Both models also use identical definitions for the optical properties of a single canopy layer as a function of leaf and wood optical properties, leaf angle distribution, canopy clumping, solar zenith angle, and leaf, wood, and crown area indices. Both models represent direct (a.k.a. “beam”) radiation as an exponentially decaying process, and solve for the diffuse (a.k.a. “hemispherical”, “isotropic”) radiation at each layer using a linear matrix equation. Where the models differ is in the terms of this matrix equation. The two-stream model (default) uses the two-stream approximation originally derived for atmospheric radiative transfer (Meador and Weaver 1980) and later adapted to vegetation canopies (Dickinson 1983, Sellers 1985). Meanwhile, the multiple-scatter model was derived from first principles by Zhao and Qualls (2005) specifically for vegetation canopies, specifically to address known biases and limitations of the two-stream approach (e.g. Wang 2003). Both of these models are described in detail in the Supplementary Information.

The crown area submodel in ED-2.2 determines the nature of competition for light between cohorts. In the default configuration (“infinite crown area”, or “complete shading”), the leaf area of a cohort is distributed across the entire horizontal area of a patch. This means that taller cohorts always receive more incoming radiation than shorter cohorts, even when the height difference is small. This has been shown to excessively suppress competition from sub-dominant individuals and result in unrealistically homogeneous canopies (Fisher et al. 2015). In the alternate configuration (“finite crown area”, or “partial shading”), canopies take up only a fraction of the available horizontal area, meaning that multiple cohorts of similar height can receive the same level of light. The horizontal area of crowns is determined by allometric equations from Dietze and Clark (2008).

The third submodel we evaluated was trait plasticity. In the default configuration, all cohorts of a given plant functional type will have the same parameters, regardless of environmental conditions. This ignores the globally-documented intraspecific trait variability as a function of light level (Niinemets 2010, Keenan and Niinemets 2016). In the alternate configuration, as light level decreases (trees become more shaded), specific leaf area increases and Vcmax decreases, following empirical relationships from the tropics (Lloyd et al. 2010).

These sub-models are summarized in Table 1.

Table 1: ED-2.2 submodel descriptions

Name Description Color
Crown model
closed (default) Cohort crowns take up entire patch area. Competition for light based only on height. Light
finite Cohort crown area is proportional to DBH according to PFT-specific allometry. Dark
Radiative transfer model
two-stream (default) Two-stream approximation (Sellers 1985, Oleson et al. 2013) Primary (red, blue)
multi-scatter Multiple-scatter approximation, following (Zhao and Qualls 2005) Secondary (green, orange)
Trait plasticity
static (default) SLA and Vc,max are constant Cool (blue, green)
plastic SLA increases, and Vc,max decreases, with light level Warm (red, orange)

Parameter descriptions

The full list of parameters in ED-2.2 is large, with over 100 parameters per plant functional type. A full sensitivity and uncertainty analysis across all of these parameters is outside the scope of this study. Instead, we selected a subset of parameters guided by (but expanding upon) previous ED2 sensitivity studies (Dietze et al. 2014, Raczka et al. 2018, Viskari et al. 2019).

Three parameters are related to leaf-level physiology: Following the enzyme-kinetic model of Farquhar et al. (1980), the rate of photosynthesis is the minimum of light-limited and enzyme-limited reactions. The former are controlled by the quantum efficiency parameter—maximum efficiency with which absorbed photosynthetically active radiation is converted to CO2. The latter are controlled by Vcmax, the maximum rate of carbon fixation by Rubisco. The water demand of photosynthesis is controlled by the stomatal slope, the sensitivity of stomatal conductance of CO2 as a function of CO2 concentration and humidity at the leaf surface (Ball et al. 1987, Leuning 1995). Three more parameters correspond to the respiration rates of leaves, roots, and “growth maintenance”.

Two more parameters control carbon allocation: One is the ratio of fine root to leaf biomass (fineroot2leaf, or q), and another is the ratio of “storage” carbon allocated to reproduction (r_fract).

Three parameters control various aspects of adult tree mortality. The overall adult mortality rate in ED-2.2 (\(M\)) is the sum of density-independent mortality from aging (\(M_I\), plants year-1), density-dependent mortality from C starvation (\(M_D\)), and mortality from cold/frost (\(M_F\)) (ED-2.2 technically also includes additional term for fire mortality, but we did not include fire in our simulations):

\[ M = M_I + M_D + M_F \]

Density-independent mortality from aging (\(M_I\)) is a prescribed, PFT-specific parameter (mort3). Density-dependent mortality from C starvation is calculated as a function of a cohort’s C balance limitation:

\[ M_D = \frac{y_1}{1 + \exp\left[y_2 \left( \frac{C_k}{C_k^*} \right)\right]} \]

where \(C_k\) is the 12 month running mean C balance of the cohort, \(C_k^*\) is the running mean ideal C balance if there was no light or water limitation, and \(y_1\) (mort1; plants year-1) and \(y_2\) (mort2; unitless) are PFT-specific parameters. Seedling mortality rate is prescribed as its own PFT-specific parameter.

Several parameters are related to canopy structure and radiative transfer. Specific leaf area (SLA) is used to convert leaf biomass to leaf area index, which in turn is used in a variety of calculations related to canopy radiative transfer and micrometeorology. Canopy clumping factor describes how evenly leaf area is distributed in horizontal space (1 being perfectly evenly; 0 being a “black hole” where all leaves are concentrated in a single point); and leaf orientation factor describes the average distribution of leaf angles (-1 being perfectly vertical, 1 being perfectly horizontal, and 0 being random). Four parameters control leaf optical properties, namely the fractions of light reflected or transmitted in visible and near-infrared wavelengths (leaf_(reflect|trans)_(vis|nir)). Details of how these parameters influence canopy radiative transfer are described in the Supplementary Information.

Other parameters:

  • Fraction of litter carbon that enters the labile (fast) soil carbon pool (f_labile)
  • Minimum plant height for reproduction (minimum_height)
  • Soil-plant water conductance, or “water availability factor” (water_conductance)
  • C:N ratio in fine roots (c2n_fineroot) and leaves (c2n_leaf)
  • Leaf turnover rate (evergreen-only; leaf_turnover_rate)
  • Proportion of dispersal that is “global” (nonlocal_dispersal)

The full list of parameters is shown in Table 2.

Table 2: Parameter descriptions. For additional information, see the ED2 model wiki.
ED Name Display name Description unit_markdown Raczka Dietze
c2n_fineroot C:N root C:N ratio in fine roots unitless (mass ratio) NA NA
c2n_leaf C:N leaf C:N ratio in leaves unitless (mass ratio) NA NA
clumping_factor Clumping Canopy clumping factor NA NA NA
f_labile F. labile Fraction of litter that goes to labile (fast) C pool unitless (0-1) Labile carbon NA
fineroot2leaf Root:leaf Ratio of fine root to leaf biomass (q) unitless (mass ratio) Root/Leaf carbon Leaf:Root
growth_resp_factor Growth resp. Fraction of daily C gain lost to growth respiration unitless (0-1) Growth respiration Growth Resp
leaf_reflect_nir Refl. (NIR) Leaf reflectance in NIR range (700-2500 nm) unitless (0-1) NA NA
leaf_reflect_vis Refl. (VIS) Leaf reflectance in visible range (400-700 nm) unitless (0-1) NA NA
leaf_respiration_rate_m2 Leaf resp. Ratio of leaf respiration to Vcmax (Rd0) ??? Leaf respiration* Leaf Resp*
leaf_trans_nir Trans. (NIR) Leaf transmittance in NIR range (700-2500 nm) unitless (0-1) NA NA
leaf_trans_vis Trans. (VIS) Leaf transmittance in visible range (400-700 nm) unitless (0-1) NA NA
leaf_turnover_rate Leaf turnover Temperature dependent rate of leaf loss (conifer only) year-1 NA NA
minimum_height Min. height Minimum height for plant reproduction m Minimum height NA
mort1 Mort. C rate Time-scale at which low-carbon balance plants die years-1 Carbon balance mortality Mortality
mort2 Mort. C bal. C balance ratio at which mortality rapidly increases unitless (mass ratio) NA NA
mort3 Mort. bg. rate Density-independent (background) mortality rate year-1 Background mortality NA
nonlocal_dispersal Dispersal Proportion of dispersal that is global unitless (0-1) NA NA
orient_factor Leaf orient. Leaf angle orientation distribution unitless (-1-1) NA NA
quantum_efficiency Quant. eff. Farquhar model parameter (TODO) mol CO2 (mol photons)-1 Quantum efficiency Quantum Eff.
r_fract Repro. C frac. Fraction of C storage to seed reproduction unitless Recruitment carbon Reproduction?
root_respiration_rate Root resp. Root respiration rate at 15 °C (root_respiration_rate/factor) μmol CO2 (kg fineroot)-1 Root respiration NA
root_turnover_rate Root turnover Temperature dependent rate of fine root loss year-1 Root turnover Root turnover
seedling_mortality Mort. seedling Proportion of seed that dies and goes to litter pool unitless (0-1) NA NA
SLA SLA Specific leaf area m2 kg-1 C Specific leaf area SLA
stomatal_slope Stomatal slope Slope between A and stomatal conductance (Leuning) unitless Stomatal sensitivity Stomatal Slope
Vcmax Vcmax Maximum rate of CO2 carboxylation at 15 °C (Vm0) μmol m-2 s-1 Vcmax Vcmax
water_conductance Water cond. Water availability factor m-2 a-1 (kg C root)-1 Soil-plant water conductance Water Cond

Plant functional types and parameterization

By default, ED-2.2 supports 17 different plant functional types, which divide plant species according to photosynthetic pathway (C3 vs. C4), growth form (grass vs. tree), leaf phenology habit (deciduous vs. evergreen), biome (e.g. temperate vs. tropical), and successional status (e.g. early, mid, late). However, we limited our simulations to the four plant functional types that have any appreciable presence at UMBS: Early, mid, and late temperate deciduous trees and pines. The species comprising these plant functional types are shown in Table 3.

Table 3: Plant functional type-species mappings

Plant functional type Species Color
Early hardwood Betula papyrifera Violet
Populus grandidentata
Populus tremuloides
Mid hardwood Quercus rubra Blue
Acer rubrum
Acer pensylvaticum
Late hardwood Acer saccharum Green
Fagus grandifolia
Pine Pinus strobus Yellow

For each plant functional type, we generated a distribution of parameter values via the Predictive Ecosystem Analyzer (PEcAn) trait-meta analysis (LeBauer et al. 2013, see also Dietze et al. 2014, Raczka et al. 2018). Prior distributions for this meta-analysis were largely adapted from previous ED2 parameter uncertainty studies (Dietze et al. 2014, Raczka et al. 2018, Viskari et al. 2019), and are shown in detail in Supplementary Information X. Species trait data for this meta-analysis came from existing records in the BETY database (www.betydb.org, LeBauer et al. 2017), as well as from publicly available records in the TRY database (www.try-db.org, Kattge et al. 2011) and, specifically for leaf optical properties, from Shiklomanov (2019 {dissertation, chapter 3}). The resulting parameter distributions are shown in Figure 1. All parameters not described in this section were set to their ED-2.2 PFT-specific default values.

**Figure 1**: Input parameter distributions from PEcAn trait meta-analysis.

Figure 1: Input parameter distributions from PEcAn trait meta-analysis.

Modeling workflow

For each factorial combination of ED-2.2 configurations (described above), we ran 100 ensemble members from 1901 to 2000. Each ensemble member was initialized from a “near-bare ground” condition: An equal number of seedlings of each plant functional type (see previous section) at the minimum resolvable size. Driving meteorological data was 6-hourly CRU-NCEP combined with an annual atmospheric CO2 record from Law-Dome ice core (Etheridge et al. 1998) and Mauna Loa observatory (Thoning et al. 1989). Soil texture was set to 92% sand, 7% silt, and 1% clay, per site-level observations in (Gough et al. 2010). The initial soil moisture profile was set to the average soil moisture profile reported in the UMBS Ameriflux ancillary data (https://ameriflux.lbl.gov/sites/siteinfo/US-UMd).

Driver uncertainty is outside the scope of this study, and initial condition uncertainty is minimized by our experimental design, which is conditioned on a specific initial condition (near bare ground).

Analysis of results

  • Overall summary:
    • Overall parameter uncertainty – Time-series plots of total annual GPP and NPP, and maximum LAI and AGB, averaged across ensemble members
    • Impact of parameter uncertainty on ecological trajectories – plot of max annual LAI by PFT
    • Both compared against observations of recently observed LAI and NPP at UMBS (e.g. FASET, FoRTE)
  • Uncertainty analysis:
    • Parameter vs. structural uncertainty – Variance in mean annual GPP from 19XX to 19YY across ensembles within model structure, and compare against
    • Breakdown of parameter uncertainty – PEcAn sensitivity and uncertainty analysis (LeBauer et al. 2013, Dietze et al. 2014)
  • Details of structural uncertainty:
    • Vertical profiles of radiation, traits?

All analyses for this work were performed using R 3.6 (see “Colophon” for full details). Code and supporting data for reproducing this analysis are publicly available on the Open Science Framework (OSF) at https://osf.io/dznuf/.

Results

Summary

**Figure 2**: ED2 plot-level predictions of gross (GPP) and net (NPP) primary productivity, aboveground biomass (AGB), total leaf area index (LAI), and Shannon diversity index (of plant functional types) by model configuration. The solid line is the ensemble mean. The dashed line is the 80% confidence interval. Black dot with error bars is the observed mean and min/max value (TODO REF).

Figure 2: ED2 plot-level predictions of gross (GPP) and net (NPP) primary productivity, aboveground biomass (AGB), total leaf area index (LAI), and Shannon diversity index (of plant functional types) by model configuration. The solid line is the ensemble mean. The dashed line is the 80% confidence interval. Black dot with error bars is the observed mean and min/max value (TODO REF).

Model projections of gross and net primary productivity, aboveground biomass, leaf area index, and functional diversity varied with both parameter value and model structure (Figure 2). All model configurations exhibited a similar trajectory of rapid growth in the first 15-20 years followed by a gradual decline to a lower, stable value around 1975 (with individual ensemble members sometimes undergoing sudden collapse in specific years). Two model configurations – finite canopy radius and two-stream radiative transfer (regardless of trait plasticity) – had much lower productivity, biomass, and diversity than the other configurations, which generally had similar ensemble mean values and spread. The ensemble means of all model configurations under-predicted observed net primary productivity and leaf area index, both at the peak of initial growth (around 1920) and at the end of the simulations, though peak predicted values were much closer to thte observations. However, the ensemble spread around the mean was large: across all configurations (except the finite canopy - two-stream radiative transfer cases), the spread in net primary productivity was from zero to nearly twice as high as the observation. Ensemble variance in leaf area index was lower on average but, in some cases (especially the combination of finite canopy and multiple-scatter radiative transfer) was over an order of magnitude. The predicted Shannon diversity index of plant functional types was, on average, higher in simulations when the canopy architecture was finite rather than closed.

**Figure 3**: Post-1975 growing season averaged NPP vs. LAI, with observation (in black) and linear regression, by model type. Each point indicates a single ensemble run with a different parameter combination. Letter labels are selected runs with the same parameters, and match the rows in Figure 4.

Figure 3: Post-1975 growing season averaged NPP vs. LAI, with observation (in black) and linear regression, by model type. Each point indicates a single ensemble run with a different parameter combination. Letter labels are selected runs with the same parameters, and match the rows in Figure 4.

Model configuration affected the emergent relationship between net primary productivity and leaf area index (Figure 3). The crown model (closed vs. finite) had the strongest effect: Compared to the simulations with “closed” crowns, simulations with finite canopy radius generally had lower net primary productivity for the same leaf area index. Although many simulations matched either the observed net primary productivity or the observed leaf area index, only a few individual runs matched both.

**Figure 4**: PFT-level predictions of leaf area index (LAI) by model configuration and parameterization. Each row of plots is a simulation with the same set of parameters across model configurations. Letters correspond to the labels in Figure 3. Note that these are predictions of total leaf area index summed across all cohorts of the same PFT, meaning that higher LAI in one PFT does not necessarily indicate dominance, size, etc.

Figure 4: PFT-level predictions of leaf area index (LAI) by model configuration and parameterization. Each row of plots is a simulation with the same set of parameters across model configurations. Letters correspond to the labels in Figure 3. Note that these are predictions of total leaf area index summed across all cohorts of the same PFT, meaning that higher LAI in one PFT does not necessarily indicate dominance, size, etc.

Model structure and parameters interacted to determine the community assembly and dynamics of each simulation, with consequences for overall ecosystem productivity that varied by model structure (Figure 4). Consider simulations A and C, the former showing competition of early hardwood with pine, and the latter completely dominated by pine: With a closed canopy and static traits, both simulations produced similar overall LAI and NPP values (Figure 3), but moving to a finite canopy dramatically increased the total LAI in C without much affecting A, while enabling trait plasticity dramatically increased the NPP of A without much affecting C (Figure 3). This example also shows the impacts of different configurations on predicted community dynamics – namely, in this specific case, enabling trait plasticity helped pine out-perform early hardwood, but only under a closed (rather than finite) canopy. Finite canopy radius similarly benefits early hardwood in its competition against late hardwood in simulations B and, to a lesser extent, D. Simulations B and D also show that the sensitivity of a given parameter combination to model structure is itself dependent on the parameters – simulation B produced similar ecological dynamics and aggregate state variables regardless of model configuration (except the extreme cases; Figure 3), while simulation D had near-zero NPP and LAI in two cases (closed, two-stream, static and finite, multi-scatter, plastic), comparable NPP and LAI to B with closed canopy and trait plasticity, and higher NPP and LAI than D with a finite canopy, multi-scatter radiative transfer, and static traits (Figure 3).

Uncertainty analysis

**Figure 5**: Comparison of parametric uncertainty within ensembles (colored bars, colored by model type) and "structural" uncertainty (variance in ensemble means; black bar) by output variable. Output variables are expressed as growing season averages for all years after 1910.

Figure 5: Comparison of parametric uncertainty within ensembles (colored bars, colored by model type) and “structural” uncertainty (variance in ensemble means; black bar) by output variable. Output variables are expressed as growing season averages for all years after 1910.

Except for the configurations with persistently low productivity, parameter uncertainty (across ensembles within each model configuration) was much larger than structural uncertainty (variance in ensemble means across configurations; Figure 4). The relative importance of parameter uncertainty varied with both model configuration and the variable in question. Allowing SLA and Vcmax to vary with light level (warm colors, as opposed to static traits – cool colors) slightly reduced variability in primary productivity, slightly increased variability in Shannon diversity, and had mixed effects on aboveground biomass and total leaf area. Enabling the finite canopy radius (dark colors; compared to closed canopy – light colors) dramatically decreased variability in primary productivity and Shannon diversity and increased variability in leaf area, and had mixed effects on aboveground biomass. reduced variance in aboveground biomass, diversity, and primary productivity. Finally, using a multiple-scatter (secondary colors) rather than a two-stream (primary colors) radiative transfer scheme increased variance in leaf area index, aboveground biomass, and diversity, but did not have a consistent effect on variance in productivity productivity.

**Figure 6**: PEcAn-like parameter sensitivity and uncertainty analysis, by model type, for 1975-2000 averaged growing season net primary productivity (left) and leaf area index (right). Elasticity (top) is the normalized sensitivity of the model to a fixed change in the parameter. Partial variance (bottom) describes the overall contribution of the parameter to model predictive uncertainty based on the combination of parameter uncertainty and model sensitivity. Each panel shows the top 10 parameter values, in terms of the absolute value of the corresponding quantity. Input parameter distributions are shown in Figure 1. Note that x-axis scales vary across panels.

Figure 6: PEcAn-like parameter sensitivity and uncertainty analysis, by model type, for 1975-2000 averaged growing season net primary productivity (left) and leaf area index (right). Elasticity (top) is the normalized sensitivity of the model to a fixed change in the parameter. Partial variance (bottom) describes the overall contribution of the parameter to model predictive uncertainty based on the combination of parameter uncertainty and model sensitivity. Each panel shows the top 10 parameter values, in terms of the absolute value of the corresponding quantity. Input parameter distributions are shown in Figure 1. Note that x-axis scales vary across panels.

Model sensitivity to specific parameters (“elasticity” sensu LeBauer et al. 2013) varied depending on the model configuration and the target output variable (Figure 5). The parameter with the single highest elasticity to NPP and LAI across all model configurations was specific leaf area (SLA), though the PFT whose SLA was most important varied with model configuration. Other frequently (but not always) important parameters were the leaf C:N ratio (c2n_leaf), ratio of fine root to leaf biomass (fineroot2leaf), leaf and root respiration rates, and leaf optical properties (leaf_reflect_vis/nir).

The parameters that contributed most to overall model predictive uncertainty (“partial variance”) differed from parameters to which the model was most sensitive (Figure 5), largely driven by the fact that many of the parameters with the highest sensitivity (e.g. SLA, Vcmax) were also well constrained by data (Figure 1). Leaf optical properties of at least one PFT were among the top 10 most important parameters across all configurations, and in multiple cases, were the most important. Other often important parameters were plant-soil water conductance (water_conductance), parameters related to adult mortality (mort1,2,3), and growth respiration. Again, patterns in partial variance were largely idiosyncratic across model configurations.

Discussion

The overarching objective of this study was to investigate how model structure and parameters interact in multidecadal simulations of forest ecosystem productivity and community assembly.

ED2.2 vs. reality

Our first question was: To what extent can the model, under any configuration, reproduce the observed reality? Ensemble mean trajectory matches traditional understanding of succession—rapid growth to peak productivity, followed by gradual decline (Odum 1969)—but happens too quickly (Figure 2); real ecosystems sustain peak productivity for longer (Gough et al. 2016). This idea that ED2.2’s main issue is that it is moving too quickly is supported by the fact that ED2.2 generally underestimates productivity after 100 years, but model predictions at peak productivity are close to reality. One reason ED2.2 is moving too quickly through succession is missing mechanisms for sustaining productivity through early succession. In particular, recent evidence points to the important role of low-severity, high-frequency disturbances in sustaining ecosystem productivity (Gough et al. 2016). ED2.2 has a mortality scheme, and indeed, parameters related to mortality were often among the top 10 most important parameters for predictive uncertainty (Figure 6). However, our simulations did not include other sources of disturbance, such as partial harvests, extreme weather, and damage from insects and pathogens, which are common to this region (Amiro et al. 2010) and which may contribute to sustaining productivity and dominance of early-successional forests (Gough et al. 2016). Future work will: (1) Calibrate ED2.2 disturbance response to experimental observations (FoRTE), and (2) run long-term simulations with prescribed disturbance frequencies and see how they affect ecological dynamics.

  • Figure 3:
    • For the same LAI, closed canopy had higher productivity than finite. Why?
      • Not zero-sum; finite canopy hurts dominant trees but doesn’t improve subdominant species that much
    • Importance of using multiple observational constraints. Ecological equifinality is a real problem.

Drivers of uncertainty

Parameter uncertainty was the dominant driver of uncertainty in our analysis. Different parameter combinations led to huge ranges in ensemble predictions of ensemble variables – NPP ranged from 0 to almost 2x the observed, and LAI was up to 3x observed (though more tightly clustered around zero) (Figure 2). Moreover, different parameter combinations produced qualitatively different ecological narratives (Figure 4). Parameter uncertainty has implications beyond C cycling—differences in species composition and relative abundance have implications for forest management and wildlife conservation, among others.

Parameter uncertainty is not synonymous with sensitivity. Many of the parameters to which model outputs are most sensitive are already well constrained by data and do not contribute as much to predictive uncertainty (Figure 6). Although parameter uncertainty is still large, this means that targeted data collection works and that we are moving in the right direction! This speaks to the value of looking at not just sensitivity, but partial variance (LeBauer et al. 2013, Dietze et al. 2014). (Compare sensitivity results to these studies, and also to Pandit et al. (2018)).

Although parameter uncertainty dominated, structural and parameter uncertainty interacted in important ways. The relative position of different parameter combinations in multivariate space (Figure 3), the total amount of parameter uncertainty (Figure 5), and the relative importance of different parameters to predictive uncertainty (Figure 6) all differed across model structures (Figure 3). This is despite the fact that our model structure changes were relatively modest compared to differences between different models. Collectively, this means that sensitivity and uncertainty analyses are not necessarily generalizable across models.

Conclusions

Acknowledgements

This project funded by NSF grant. Cyberinfrastructure provided by Pacific Northwest National Laboratory (PNNL). Data from University of Michigan Biological Station (UMBS). Data from TRY (TODO: Specific TRY statement).

pagebreak

References

Amiro, B. D., A. G. Barr, J. G. Barr, T. A. Black, R. Bracho, M. Brown, J. Chen, K. L. Clark, K. J. Davis, A. R. Desai, and et al. 2010. Ecosystem carbon dioxide fluxes after disturbance in forests of North America. Journal of Geophysical Research 115.

Ball, J. T., I. E. Woodrow, and J. A. Berry. 1987. A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions. Progress in Photosynthesis Research:221–224.

Bergen, K. M., and I. Dronova. 2007. Observing succession on aspen-dominated landscapes using a remote sensing-ecosystem approach. Landscape Ecology 22:1395–1410.

Birdsey, R., K. Pregitzer, and A. Lucier. 2006. Forest carbon management in the united states. Journal of Environment Quality 35:1461.

De Kauwe, M. G., B. E. Medlyn, S. Zaehle, A. P. Walker, M. C. Dietze, T. Hickler, A. K. Jain, Y. Luo, W. J. Parton, I. C. Prentice, and et al. 2013. Forest water use and water use efficiency at elevated CO2: A model-data intercomparison at two contrasting temperate forest FACE sites. Global Change Biology 19:1759–1779.

Dickinson, R. E. 1983. Land surface processes and climate-surface albedos and energy balance. Theory of Climate, Proceedings of a Symposium Commemorating the Two-Hundredth Anniversary of the Academy of Sciences of Lisbon:305–353.

Dietze, M. C., and J. S. Clark. 2008. Changing the gap dynamics paradigm: Vegetative regeneration control on forest response to disturbance. Ecological Monographs 78:331–347.

Dietze, M. C., S. P. Serbin, C. Davidson, A. R. Desai, X. Feng, R. Kelly, R. Kooper, D. LeBauer, J. Mantooth, K. McHenry, and et al. 2014. A quantitative assessment of a terrestrial biosphere model’s data needs across North American biomes. Journal of Geophysical Research: Biogeosciences 119:286–300.

Etheridge, D. M., L. P. Steele, R. L. Langefelds, R. J. Francey, J.-M. Marnola, and V. I. Morgan. 1998. Historical CO2 records from the Law Dome DE08, DE08-2, and DSS ice cores. In trends: A compendium of data on global change. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA.

Farquhar, G. D., S. von Caemmerer, and J. A. Berry. 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149:78–90.

Fisher, R. A., C. D. Koven, W. R. L. Anderegg, B. O. Christoffersen, M. C. Dietze, C. E. Farrior, J. A. Holm, G. C. Hurtt, R. G. Knox, P. J. Lawrence, J. W. Lichstein, M. Longo, A. M. Matheny, D. Medvigy, H. C. Muller-Landau, T. L. Powell, S. P. Serbin, H. Sato, J. K. Shuman, B. Smith, A. T. Trugman, T. Viskari, H. Verbeeck, E. Weng, C. Xu, X. Xu, T. Zhang, and P. R. Moorcroft. 2017. Vegetation demographics in earth system models: A review of progress and priorities. Global Change Biology 24:35–54.

Fisher, R. A., S. Muszala, M. Verteinstein, P. Lawrence, C. Xu, N. G. McDowell, R. G. Knox, C. Koven, J. Holm, B. M. Rogers, and et al. 2015. Taking off the training wheels: The properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED). Geoscientific Model Development 8:3593–3619.

Fisher, R., N. McDowell, D. Purves, P. Moorcroft, S. Sitch, P. Cox, C. Huntingford, P. Meir, and F. Ian Woodward. 2010. Assessing uncertainties in a second-generation dynamic vegetation model caused by ecological scale limitations. New Phytologist 187:666–681.

Franks, P. J., G. B. Bonan, J. A. Berry, D. L. Lombardozzi, N. M. Holbrook, N. Herold, and K. W. Oleson. 2018. Comparing optimal and empirical stomatal conductance models for application in Earth system models. Global Change Biology 24:5708–5723.

Friedlingstein, P., P. Cox, R. Betts, L. Bopp, W. von Bloh, V. Brovkin, P. Cadule, S. Doney, M. Eby, I. Fung, and et al. 2006. Climate-carbon cycle feedback analysis: Results from the C4MIP model intercomparison. Journal of Climate 19:3337–3353.

Friedlingstein, P., M. Meinshausen, V. K. Arora, C. D. Jones, A. Anav, S. K. Liddicoat, and R. Knutti. 2014. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. Journal of Climate 27:511–526.

Gough, C. M., P. S. Curtis, B. S. Hardiman, C. M. Scheuermann, and B. Bond-Lamberty. 2016. Disturbance, complexity, and succession of net ecosystem production in North America’s temperate deciduous forests. Ecosphere 7:e01375.

Gough, C. M., C. S. Vogel, B. Hardiman, and P. S. Curtis. 2010. Wood net primary production resilience in an unmanaged forest transitioning from early to middle succession. Forest Ecology and Management 260:36–41.

Gustafson, E. J., A. M. G. De Bruijn, B. R. Miranda, and B. R. Sturtevant. 2016. Implications of mechanistic modeling of drought effects on growth and competition in forest landscape models. Ecosphere 7:e01253.

Hawkes, C. 2000. Woody plant mortality algorithms: Description, problems and progress. Ecological Modelling 126:225–248.

Hurtt, G. C., R. Dubayah, J. Drake, P. R. Moorcroft, S. W. Pacala, J. B. Blair, and M. G. Fearon. 2004. Beyond potential vegetation: Combining lidar data and a height-structured model for carbon studies. Ecological Applications 14:873–883.

Kattge, J., S. Díaz, S. Lavorel, I. C. Prentice, P. Leadley, G. Bönisch, E. Garnier, M. Westoby, P. B. Reich, I. J. Wright, J. H. C. Cornelissen, C. Violle, S. P. Harrison, P. M. Van Bodegom, M. Reichstein, B. J. Enquist, N. A. Soudzilovskaia, D. D. Ackerly, M. Anand, O. Atkin, M. Bahn, T. R. Baker, D. Baldocchi, R. Bekker, C. C. Blanco, B. Blonder, W. J. Bond, R. Bradstock, D. E. Bunker, F. Casanoves, J. Cavender-Bares, J. Q. Chambers, F. S. Chapin III, J. Chave, D. Coomes, W. K. Cornwell, J. M. Craine, B. H. Dobrin, L. Duarte, W. Durka, J. Elser, G. Esser, M. Estiarte, W. F. Fagan, J. Fang, F. Fernández-Méndez, A. Fidelis, B. Finegan, O. Flores, H. Ford, D. Frank, G. T. Freschet, N. M. Fyllas, R. V. Gallagher, W. A. Green, A. G. Gutierrez, T. Hickler, S. I. Higgins, J. G. Hodgson, A. Jalili, S. Jansen, C. A. Joly, A. J. Kerkhoff, D. Kirkup, K. Kitajima, M. Kleyer, S. Klotz, J. M. H. Knops, K. Kramer, I. Kühn, H. Kurokawa, D. Laughlin, T. D. Lee, M. Leishman, F. Lens, T. Lenz, S. L. Lewis, J. Lloyd, J. Llusià, F. Louault, S. Ma, M. D. Mahecha, P. Manning, T. Massad, B. E. Medlyn, J. Messier, A. T. Moles, S. C. Müller, K. Nadrowski, S. Naeem, Ü. Niinemets, S. Nöllert, A. Nüske, R. Ogaya, J. Oleksyn, V. G. Onipchenko, Y. Onoda, J. Ordoñez, G. Overbeck, W. A. Ozinga, S. Patiño, S. Paula, J. G. Pausas, J. Peñuelas, O. L. Phillips, V. Pillar, H. Poorter, L. Poorter, P. Poschlod, A. Prinzing, R. Proulx, A. Rammig, S. Reinsch, B. Reu, L. Sack, B. Salgado-Negret, J. Sardans, S. Shiodera, B. Shipley, A. Siefert, E. Sosinski, J.-F. Soussana, E. Swaine, N. Swenson, K. Thompson, P. Thornton, M. Waldram, E. Weiher, M. White, S. White, S. J. Wright, B. Yguel, S. Zaehle, A. E. Zanne, and C. Wirth. 2011. TRY – a global database of plant traits. Global Change Biology 17:2905–2935.

Keenan, T. F., and Ü. Niinemets. 2016. Global leaf trait estimates biased due to plasticity in the shade. Nature Plants 3.

LeBauer, D., R. Kooper, P. Mulrooney, S. Rohde, D. Wang, S. P. Long, and M. C. Dietze. 2017. BETYdb: A yield, trait, and ecosystem service database applied to second-generation bioenergy feedstock production. GCB Bioenergy 10:61–71.

LeBauer, D. S., D. Wang, K. T. Richter, C. C. Davidson, and M. C. Dietze. 2013. Facilitating feedbacks between field measurements and ecosystem models. Ecological Monographs 83:133–154.

Leuning, R. 1995. A critical appraisal of a combined stomatal-photosynthesis model for C3 plants. Plant, Cell and Environment 18:339–355.

Liénard, J., J. Harrison, and N. Strigul. 2016. US forest response to projected climate-related stress: Atoleranceperspective. Global Change Biology 22:2875–2886.

Lloyd, J., S. Patiño, R. Q. Paiva, G. B. Nardoto, C. A. Quesada, A. J. B. Santos, T. R. Baker, W. A. Brand, I. Hilke, H. Gielmann, and et al. 2010. Optimisation of photosynthetic carbon gain and within-canopy gradients of associated foliar traits for Amazon forest trees. Biogeosciences 7:1833–1859.

Longo, M., R. G. Knox, N. M. Levine, A. L. S. Swann, D. M. Medvigy, M. C. Dietze, Y. Kim, K. Zhang, D. Bonal, B. Burban, and et al. 2019a. The biophysics, ecology, and biogeochemistry of functionally diverse, vertically- and horizontally-heterogeneous ecosystems: The Ecosystem Demography model, version 2.2 - part 2: Model evaluation. Geoscientific Model Development Discussions:1–34.

Longo, M., R. G. Knox, D. M. Medvigy, N. M. Levine, M. C. Dietze, Y. Kim, A. L. S. Swann, K. Zhang, C. R. Rollinson, R. L. Bras, and et al. 2019b. The biophysics, ecology, and biogeochemistry of functionally diverse, vertically- and horizontally-heterogeneous ecosystems: The Ecosystem Demography model, version 2.2 - part 1: Model description. Geoscientific Model Development Discussions:1–53.

Lovenduski, N. S., and G. B. Bonan. 2017. Reducing uncertainty in projections of terrestrial carbon uptake. Environmental Research Letters 12:044020.

Lovett, G. M., M. Weiss, A. M. Liebhold, T. P. Holmes, B. Leung, K. F. Lambert, D. A. Orwig, F. T. Campbell, J. Rosenthal, D. G. McCullough, and et al. 2016. Nonnative forest insects and pathogens in the United States: Impacts and policy options. Ecological Applications 26:1437–1455.

Meador, W. E., and W. R. Weaver. 1980. Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement. Journal of the Atmospheric Sciences 37:630–643.

Medlyn, B. E., S. Zaehle, M. G. De Kauwe, A. P. Walker, M. C. Dietze, P. J. Hanson, T. Hickler, A. K. Jain, Y. Luo, W. Parton, and et al. 2015. Using ecosystem experiments to improve vegetation models. Nature Climate Change 5:528–534.

Medvigy, D., and P. R. Moorcroft. 2011. Predicting ecosystem dynamics at regional scales: An evaluation of a terrestrial biosphere model for the forests of northeastern North America. Philosophical Transactions of the Royal Society B: Biological Sciences 367:222–235.

Medvigy, D., S. C. Wofsy, J. W. Munger, D. Y. Hollinger, and P. R. Moorcroft. 2009. Mechanistic scaling of ecosystem function and dynamics in space and time: Ecosystem demography model version 2. Journal of Geophysical Research 114.

Miller, A. D., M. C. Dietze, E. H. DeLucia, and K. J. Anderson-Teixeira. 2015. Alteration of forest succession and carbon cycling under elevated CO2. Global Change Biology 22:351–363.

Mladenoff, D. J. 2004. LANDIS and forest landscape models. Ecological Modelling 180:7–19.

Moorcroft, P. R., G. C. Hurtt, and S. W. Pacala. 2001. A method for scaling vegetation dynamics: The ecosystem demography model (ed). Ecological Monographs 71:557–586.

Niinemets, Ü. 2010. A review of light interception in plant stands from leaf to canopy in different plant functional types and in species with varying shade tolerance. Ecological Research 25:693–714.

Odum, E. P. 1969. The strategy of ecosystem development. Science 164:262–270.

Oleson, K. W., D. M. Lawrence, G. B. Bonan, B. Drewniak, M. Huang, C. D. Koven, S. Levis, F. Li, W. J. Riley, Z. M. Subin, S. C. Swenson, P. Thornton, A. Bozbiyik, R. Fisher, C. L. Heald, E. Kluzek, J.-F. Lamarque, P. J. Lawrence, R. Leung, W. Lipscom, S. Muszala, D. M. Ricciuto, W. Sacks, Y. Sun, J. Tang, and Y. Zong-Liang. 2013. Technical description of version 4.5 of the Community Land Model (CLM). NCAR Earth System Laboratory Climate; Global Dynamics Division.

Pandit, K., H. Dashti, N. F. Glenn, A. N. Flores, K. C. Maguire, D. J. Shinneman, G. N. Flerchinger, and A. W. Fellows. 2018. Optimizing shrub parameters to estimate gross primary production of the sagebrush ecosystem using the Ecosystem Demography (EDv2.2) model. Geoscientific Model Development Discussions:1–23.

Purves, D., and S. Pacala. 2008. Predictive models of forest dynamics. Science 320:1452–1453.

Raczka, B., M. C. Dietze, S. P. Serbin, and K. J. Davis. 2018. What limits predictive certainty of long‐term carbon uptake? Journal of Geophysical Research: Biogeosciences 123:3570–3588.

Rogers, A., B. E. Medlyn, J. S. Dukes, G. Bonan, S. von Caemmerer, M. C. Dietze, J. Kattge, A. D. B. Leakey, L. M. Mercado, Ü. Niinemets, and et al. 2016. A roadmap for improving the representation of photosynthesis in earth system models. New Phytologist 213:22–42.

Rollinson, C. R., Y. Liu, A. Raiho, D. J. P. Moore, J. McLachlan, D. A. Bishop, A. Dye, J. H. Matthes, A. Hessl, T. Hickler, and et al. 2017. Emergent climate and CO2 sensitivities of net primary productivity in ecosystem models do not agree with empirical data in temperate forests of eastern North America. Global Change Biology 23:2755–2767.

Scheller, R. M., J. B. Domingo, B. R. Sturtevant, J. S. Williams, A. Rudy, E. J. Gustafson, and D. J. Mladenoff. 2007. Design, development, and application of LANDIS-II, a spatial landscape simulation model with flexible temporal and spatial resolution. Ecological Modelling 201:409–419.

Sellers, P. J. 1985. Canopy reflectance, photosynthesis and transpiration. International Journal of Remote Sensing 6:1335–1372.

Shifley, S. R., F. X. Aguilar, N. Song, S. I. Stewart, D. J. Nowak, D. D. Gormanson, W. K. Moser, S. Wormstead, and E. J. Greenfield. 2012. Forests of the Northern United States. U.S. Department of Agriculture, Forest Service, Northern Research Station.

Shifley, S. R., and W. K. Moser. 2016. Future forests of the northern United States. U.S. Department of Agriculture, Forest Service, Northern Research Station.

Shividenko, A., C. V. Barber, R. Persson, P. Gonzalez, R. Hassan, P. Lakyda, I. McCallum, S. Nilsson, J. Pulhin, B. van Rosenburg, and B. Scholes. 2005. Forest and woodland ecosystems. in R. Hassan, R. Scholes, and N. Ash, editors. Millennium ecosystem assessment: Ecosystems and human well-being: Current state and trends, volume 1. Island Press.

Shugart, H. H., G. P. Asner, R. Fischer, A. Huth, N. Knapp, T. Le Toan, and J. K. Shuman. 2015. Computer and remote-sensing infrastructure to enhance large-scale testing of individual-based forest models. Frontiers in Ecology and the Environment 13:503–511.

Stearns, F., and G. E. Likens. 2002. One hundred years of recovery of a pine forest in northern Wisconsin. The American Midland Naturalist 148:2–19.

Sulman, B. N., J. A. M. Moore, R. Abramoff, C. Averill, S. Kivlin, K. Georgiou, B. Sridhar, M. D. Hartman, G. Wang, W. R. Wieder, and et al. 2018. Multiple models and experiments underscore large uncertainty in soil carbon dynamics. Biogeochemistry 141:109–123.

Swann, A. L. S., M. Longo, R. G. Knox, E. Lee, and P. R. Moorcroft. 2015. Future deforestation in the Amazon and consequences for South American climate. Agricultural and Forest Meteorology 214-215:12–24.

Swanston, C., L. A. Brandt, M. K. Janowiak, S. D. Handler, P. Butler-Leopold, L. Iverson, F. R. Thompson III, T. A. Ontl, and P. D. Shannon. 2017. Vulnerability of forests of the Midwest and Northeast United States to climate change. Climatic Change 146:103–116.

Thompson, J. R., D. N. Carpenter, C. V. Cogbill, and D. R. Foster. 2013. Four centuries of change in northeastern United States forests. PLoS ONE 8:e72540.

Thoning, K. W., P. P. Tans, and W. D. Komhyr. 1989. Atmospheric carbon dioxide at Mauna Loa Observatory: 2. Analysis of the NOAA GMCC data, 1974-1985. Journal of Geophysical Research: Atmospheres 94:8549–8565.

Viskari, T., A. Shiklomanov, M. C. Dietze, and S. P. Serbin. 2019. The influence of canopy radiation parameter uncertainty on model projections of terrestrial carbon and energy cycling. PLOS ONE 14:e0216512.

Walker, A. P., S. Zaehle, B. E. Medlyn, M. G. De Kauwe, S. Asao, T. Hickler, W. Parton, D. M. Ricciuto, Y.-P. Wang, and D. W. et al. 2015. Predicting long-term carbon sequestration in response to CO2 enrichment: How and why do current ecosystem models differ? Global Biogeochemical Cycles 29:476–495.

Wang, Y. P. 2003. A comparison of three different canopy radiation models commonly used in plant modelling. Functional Plant Biology 30:143.

Wieder, W. R., M. D. Hartman, B. N. Sulman, Y.-P. Wang, C. D. Koven, and G. B. Bonan. 2017. Carbon cycle confidence and uncertainty: Exploring variation among soil biogeochemical models. Global Change Biology 24:1563–1579.

Williams, C. A., G. J. Collatz, J. Masek, and S. N. Goward. 2012. Carbon consequences of forest disturbance and recovery across the conterminous United States. Global Biogeochemical Cycles 26:n/a–n/a.

Wolter, P. T., and M. A. White. 2002. Recent forest cover type transitions and landscape structural changes in northeast Minnesota, USA. Landscape Ecology 17:133–155.

Zaehle, S., B. E. Medlyn, M. G. De Kauwe, A. P. Walker, M. C. Dietze, T. Hickler, Y. Luo, Y.-P. Wang, B. El-Masri, P. Thornton, and et al. 2014. Evaluation of 11 terrestrial carbon-nitrogen cycle models against observations from two temperate free-air CO2 enrichment studies. New Phytologist 202:803–822.

Zhao, W., and R. J. Qualls. 2005. A multiple-layer canopy scattering model to simulate shortwave radiation distribution within a homogeneous plant canopy. Water Resources Research 41.

pagebreak

Colophon

This report was generated on 2019-09-12 13:37:40 using the following computational environment and dependencies:

#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.6.1 (2019-07-05)
#>  os       macOS High Sierra 10.13.6   
#>  system   x86_64, darwin17.7.0        
#>  ui       unknown                     
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2019-09-12                  
#> 
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package       * version    date       lib source                        
#>  assertthat      0.2.1      2019-03-21 [1] CRAN (R 3.6.1)                
#>  backports       1.1.4      2019-04-10 [1] CRAN (R 3.6.1)                
#>  base64url       1.4        2018-05-14 [1] CRAN (R 3.6.1)                
#>  callr           3.3.1      2019-07-18 [1] CRAN (R 3.6.1)                
#>  cli             1.1.0      2019-03-19 [1] CRAN (R 3.6.1)                
#>  codetools       0.2-16     2018-12-24 [3] CRAN (R 3.6.1)                
#>  colorspace      1.4-1      2019-03-18 [1] CRAN (R 3.6.1)                
#>  cowplot       * 1.0.0      2019-07-11 [1] CRAN (R 3.6.1)                
#>  crayon          1.3.4      2017-09-16 [1] CRAN (R 3.6.1)                
#>  data.table    * 1.12.2     2019-04-07 [1] CRAN (R 3.6.1)                
#>  desc            1.2.0      2018-05-01 [1] CRAN (R 3.6.1)                
#>  devtools        2.2.0      2019-09-07 [1] CRAN (R 3.6.1)                
#>  digest          0.6.20     2019-07-04 [1] CRAN (R 3.6.1)                
#>  dplyr         * 0.8.3      2019-07-04 [1] CRAN (R 3.6.1)                
#>  drake         * 7.6.1      2019-08-19 [1] CRAN (R 3.6.1)                
#>  DT              0.8        2019-08-07 [1] CRAN (R 3.6.1)                
#>  ellipsis        0.2.0.1    2019-07-02 [1] CRAN (R 3.6.1)                
#>  evaluate        0.14       2019-05-28 [1] CRAN (R 3.6.1)                
#>  filelock        1.0.2      2018-10-05 [1] CRAN (R 3.6.1)                
#>  forcats       * 0.4.0      2019-02-17 [1] CRAN (R 3.6.1)                
#>  fortebaseline * 0.0.0.9000 2019-08-09 [1] local                         
#>  fs            * 1.3.1      2019-05-06 [1] CRAN (R 3.6.1)                
#>  fst           * 0.9.0      2019-04-09 [1] CRAN (R 3.6.1)                
#>  furrr         * 0.1.0      2018-05-16 [1] CRAN (R 3.6.1)                
#>  future        * 1.14.0     2019-07-02 [1] CRAN (R 3.6.1)                
#>  future.callr  * 0.4.0      2019-01-07 [1] CRAN (R 3.6.1)                
#>  ggplot2       * 3.2.1      2019-08-10 [1] CRAN (R 3.6.1)                
#>  ggrepel       * 0.8.1      2019-05-07 [1] CRAN (R 3.6.1)                
#>  globals         0.12.4     2018-10-11 [1] CRAN (R 3.6.1)                
#>  glue            1.3.1      2019-03-12 [1] CRAN (R 3.6.1)                
#>  gtable          0.3.0      2019-03-25 [1] CRAN (R 3.6.1)                
#>  here          * 0.1        2017-05-28 [1] CRAN (R 3.6.1)                
#>  highr           0.8        2019-03-20 [1] CRAN (R 3.6.1)                
#>  hms             0.5.1      2019-08-23 [1] CRAN (R 3.6.1)                
#>  htmltools       0.3.6      2017-04-28 [1] CRAN (R 3.6.1)                
#>  htmlwidgets     1.3        2018-09-30 [1] CRAN (R 3.6.1)                
#>  igraph          1.2.4.1    2019-04-22 [1] CRAN (R 3.6.1)                
#>  knitr         * 1.24       2019-08-08 [1] CRAN (R 3.6.1)                
#>  lazyeval        0.2.2      2019-03-15 [1] CRAN (R 3.6.1)                
#>  lifecycle       0.1.0      2019-08-01 [1] CRAN (R 3.6.1)                
#>  listenv         0.7.0      2018-01-21 [1] CRAN (R 3.6.1)                
#>  lubridate     * 1.7.4      2018-04-11 [1] CRAN (R 3.6.1)                
#>  magrittr        1.5        2014-11-22 [1] CRAN (R 3.6.1)                
#>  memoise         1.1.0      2017-04-21 [1] CRAN (R 3.6.1)                
#>  munsell         0.5.0      2018-06-12 [1] CRAN (R 3.6.1)                
#>  pillar          1.4.2      2019-06-29 [1] CRAN (R 3.6.1)                
#>  pkgbuild        1.0.5      2019-08-26 [1] CRAN (R 3.6.1)                
#>  pkgconfig       2.0.2      2018-08-16 [1] CRAN (R 3.6.1)                
#>  pkgload         1.0.2      2018-10-29 [1] CRAN (R 3.6.1)                
#>  png             0.1-7      2013-12-03 [1] CRAN (R 3.6.1)                
#>  prettyunits     1.0.2      2015-07-13 [1] CRAN (R 3.6.1)                
#>  processx        3.4.1      2019-07-18 [1] CRAN (R 3.6.1)                
#>  ps              1.3.0      2018-12-21 [1] CRAN (R 3.6.1)                
#>  purrr         * 0.3.2      2019-03-15 [1] CRAN (R 3.6.1)                
#>  R6              2.4.0      2019-02-14 [1] CRAN (R 3.6.1)                
#>  Rcpp            1.0.2      2019-07-25 [1] CRAN (R 3.6.1)                
#>  readr         * 1.3.1      2018-12-21 [1] CRAN (R 3.6.1)                
#>  remotes         2.1.0      2019-06-24 [1] CRAN (R 3.6.1)                
#>  rlang           0.4.0      2019-06-25 [1] CRAN (R 3.6.1)                
#>  rmarkdown       1.15       2019-08-21 [1] CRAN (R 3.6.1)                
#>  rprojroot       1.3-2      2018-01-03 [1] CRAN (R 3.6.1)                
#>  scales          1.0.0      2018-08-09 [1] CRAN (R 3.6.1)                
#>  sessioninfo     1.1.1      2018-11-05 [1] CRAN (R 3.6.1)                
#>  storr           1.2.1      2018-10-18 [1] CRAN (R 3.6.1)                
#>  stringi         1.4.3      2019-03-12 [1] CRAN (R 3.6.1)                
#>  stringr         1.4.0      2019-02-10 [1] CRAN (R 3.6.1)                
#>  testthat        2.2.1      2019-07-25 [1] CRAN (R 3.6.1)                
#>  tibble          2.1.3      2019-06-06 [1] CRAN (R 3.6.1)                
#>  tidyr         * 1.0.0      2019-09-11 [1] CRAN (R 3.6.1)                
#>  tidyselect      0.2.5      2018-10-11 [1] CRAN (R 3.6.1)                
#>  txtq            0.1.5      2019-08-19 [1] CRAN (R 3.6.1)                
#>  usethis         1.5.1.9000 2019-09-12 [1] Github (r-lib/usethis@a2342b8)
#>  vctrs           0.2.0      2019-07-05 [1] CRAN (R 3.6.1)                
#>  withr           2.1.2      2018-03-15 [1] CRAN (R 3.6.1)                
#>  xfun            0.9        2019-08-21 [1] CRAN (R 3.6.1)                
#>  yaml            2.2.0      2018-07-25 [1] CRAN (R 3.6.1)                
#>  zeallot         0.1.0      2018-01-28 [1] CRAN (R 3.6.1)                
#> 
#> [1] /Users/shik544/R
#> [2] /usr/local/lib/R/3.6/site-library
#> [3] /usr/local/Cellar/r/3.6.1_1/lib/R/library

The current Git commit details are:

#> Local:    master /Users/shik544/Projects/forte_project/fortebaseline
#> Remote:   master @ origin (git@github.com:ashiklom/fortebaseline.git)
#> Head:     [fe7c1ec] 2019-09-12: Big-picture outline of 2nd half of discussion

View this project on GitHub in repository ashiklom/fortebaseline (this manuscript file is in analysis/paper/paper.Rmd) and in Open Science Framework (OSF) at https://osf.io/dznuf/ .